Finite temperature hydrodynamic modes of trapped quantum gases 



Andras Csordas 

Research Group for Statistical Physics of the Hungarian Academy of Sciences 
Pdzmdny Peter Setdny 1/ A 
1111 Budapest, Hungary 

Robert Graham 
Fachbereich Physik, Universitdt-GH, 
45117 Essen, Germany 

The hydrodynamic equations of an ideal fluid formed by a dilute quantum gas in a parabolic 
trapping potential are studied analytically and numerically. Due to the appearance of internal 
modes in the fluid stratified by the trapping potential, the spectrum of low-lying modes is found to 
be dense in the high-temperature limit, with an infinitely degenerate set of zero- frequency modes. 
The spectrum for Bose-fiuids and Fermi-fiuids is obtained and discussed. 



I. INTRODUCTION 
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I The successful trapping of dilute Bose- and Fermi-gases in magnetic traps and their subsequent cooling to temper- 
■ atures below quantum-degeneracy has made the study of their hydrodynamics a subject of high current interest. 
The basic hydrodynamic equations of the fluid formed by such gases in local thermodynamic equilibrium are well 
^ , known. In the limit where the fluid can be considered ideal they are the continuity equation for the mass-density p 
10 ' and velocity field u 

i> : 

the Euler equation for the velocity field 

O 

^ + u-yu^--VP + fix), (1.2) 



with the external force per unit mass / = —^'^V{x) and the pressure P related to the internal energy density e by 



P = §£. It satisfies 



8 : ap ^ „ 5 



f{uP)^pu-f (1.3) 



The thermodynamic equilibrium distributions of the density po{x) and pressure Po{x) with VPq = pof are given by 



?H ' the ideal quantum-gas expressions at constant temperature (3 = l/ksT and chemical potential 



Po{x)=m / -Jl^hi^^^) (1-4) 

/ 4^E^h{k.Sl (1.5) 
3 J (ztt)'^ 2m 



with the single-particle distribution 



/t ^) 



(1.6) 



of the Bose-Einstein (upper sign) or Fermi-Dirac (lower sign) form. 

The derivation of these equations from the Boltzmann-equation is well-known, see jj]. In recent years a number 
of papers have already been devoted to the study of solutions of these equations linearized around the equilibrium 
state in parabolic traps. Griffin, Wu and Stringari ||] derived a closed equation for the velocity fluctuations and 
gave explicit solutions for surface waves of a Bose gas in an isotropic trap and, for a classical gas, also in the axially 
symmetric anisotropic trap. In the latter case they gave also solutions for modes corresponding to irrotational flow. A 
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further studies of the hydrodynamic regime of a trapped bose-gases was presented in [|| . Fermi-gases were considered 
by Bruun and Clark [[7|. Besides considering the low temperature limit for the degenerate Fermi gas, these authors 
gave an analytical solution for the mode spectrum in an isotropic trap in the high-temperature limit and identified one 
branch of the dispersion relation as 'internal waves' driven by the inhomogeneous trap potential. This is a point which 
we intend to examine further in the present paper. Amoruso ct al. |^ also derived special solutions to the linearized 
hydrodynamic equations for the low-temperature limit of the degenerate Fermi gas. In a recent paper |^ we have 
also studied this low-temperature regime for Fermi gases and gave solutions for the completely anisotropic parabolic 
trap. The present paper will therefore concentrate on temperatures of the order of the degeneracy temperature or 
above. In a number of papers effects beyond the scope of eqs.(l.l) - (1.6) were also considered. Vichi and Stringari 
10 considered the effects of mean- fields due to interactions on the collective oscillations of Fermi gases in a trap, 
while Pethick, Smith and collaborators ||l^, |13|, Vichi and Guery-Odelin et al. discussed the coUisional 
damping of collective modes in Bose gases and Fermi gases respectively. In [|ll],0 a simple interpolation formula was 
proposed between the mode- frequencies lOc in the coUisionless regime and the hydrodynamic regime, LOh^ of the form 

(1.7) 

1 — IU}T 



where 

- < pow > (1.8) 



_i 87ra2 



is the mean collision rate. This description was further examined in ref. [[15[ . Eq.(L7) is based on general considerations 
of non-equilibrium thermodynamics | [T^ . Damping of the hydrodynamics in Bose gases was also studied by Grifhn 
and collaborators in a series of papers, see where further references can be found. 

In the present paper we will not be concerned with damping effects. Instead our goal in the present paper is to study 
further the collision-dominated dissipation-less hydrodynamic regime in harmonic traps with arbitrary anisotropy. 



We do this by giving a systematic treatment of the linearized hydrodynamic equations based on eqs.(l.l) - (1.6) 
applicable (within our basic assumptions) in the whole temperature range from the high temperature domain, where 

the Boltzmann limit /zp ~ exp(— /3 {^^^^ + l^(a') ~ /^^) applies, to the regime close to the degeneracy temperature 
for bosons and down to nearly vanishing temperature for fermions. We shall discuss a class of exact solutions of 
the dissipation-less equations applicable to the whole temperature-domain covered by the theory, generalizing results 
obtained previously for traps with axial symmetry. It can be shown that in the high-temperature limit of a classical 
Boltzmann gas the linearized hydrodynamic equations in a completely anisotropic trap are integrable and separable in 
elliptic coordinates, just as their low temperature counterparts However, at lower temperatures where effects of 

quantum statistics become important, the integrability and separability are lost, which manifests itself e.g. in effects 
of avoided level crossings. 

Of special interest in the present paper, besides the common sound modes, will be the phenomenon of 'internal 
waves', which are characteristic of fluids whose equilibrium state is stratiflcd by an external potential. Internal waves 
in trapped Fermi gases were mentioned in Q but have not yet been investigated in detail for trapped quantum-gases. 
For the discussion of internal waves in classical contexts like waves in the atmosphere see |19| . 



II. LINEARIZED HYDRODYNAMIC EQUATIONS AND HILBERT-SPACE OF THEIR SOLUTIONS 

In the present section we write down the five linearized hydrodynamic equations whose solution is the central theme 
of this paper. In previous work on these equations they were reduced to a set of three wave-equations for the velocity 
field, which we shall also write down for completeness, and some special solutions of this latter set of equations were 
given. However, the appropriate boundary conditions are hard to formulate for the velocity field and therefore it is 
not clear, so far, which function space is spanned by the solutions, and whether a scalar product can be placed on this 
function space, and if so what it is. This question is of particular practical and theoretical relevance for the present 
problem, because in general the solutions have to be constructed numerically by converting the differential operators 
to matrices using a basis and the scalar product in the solution-space. It is important to choose the correct scalar 
product because, as we shall see, the problem possesses a dense-lying discrete spectrum of low-lying states, and it is a 
priori far from clear, whether all these states are needed to span the complete space of solutions, and if not, how the 
correct states are to be distinguished. It is our aim here to devote particular attention to this open problem and to 
present an answer. The way to achieve this will be to deviate from the previous line of approach by deriving, instead 
of three coupled wave-equations for the components of the velocity field, two coupled wave-equations for the pressure 
and the density. For these we shall construct a scalar product in which the wave-operator is self-adjoint so that its 
eigenf unctions form a complete set in a well-defined Hilbert space. 
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A. Linearized hydrodynamic equations 



Let us introduce small deviations 5p, 5P of density and pressure from equilibrium 



p = po (x) + Sp (f , t),P^PQ (x) + SP {x, t) 



(2.1) 



where Pq{x) and Pq{x) solve the time-independent hydrodynamic equations with vanishing velocity- field Uq — Q, 
namely V Pq{x) = pQ{x)f{x), which defines the mechanical equilibrium condition. In principle there are many equi- 
librium profiles po{x), Po{x) satisfying this requirement. In our present context the physically relevant one is the 
thermodynamic equilibrium of maximum local entropy. The entropy-maximum is achieved by the spec ial p r ofile s 
po{x), Po{x) corresponding to a state with uniform temperature T and chemical potential p. Using eqs.(L4), ( |l.5| ) 
Pq{x), Po{x) can be written as 



rnksT 



3/2 



mkgT 
27rh^ 



3 V{x 
2 



3/2 



ksT ' 
5 V{x) - p 



kuT 



(2.2) 



where the upper (lower) sign refers to bosons (fermions) and 

A_^\ , A+=2. (2.3) 

Eqs.(^.2[), but with space-dependent T{x), p{x) apply also to the states of local thermodynamic equilibrium, in which 
the system always is in the hydrodynamic limit. They can then be taken to define two of the four fields P{x), 
p{x), T{x), p{x) in terms of the other two. In the fermionic case we need to assume the presence of two equally 
populated hyperfine sub-states, the collisions between which can then ensure the local thermodynamic equilibrium. 
The Bose-Einstein integrals F„(s,z) and Fermi-Dirac integrals F+(s,z) are defined by 



1 



T{s) 



■ dx 



(2.4) 



satisfying the familiar recursion relation. 



dz 



(2.5) 



In the present case z, and therefore also F^{s,z), is space-dependent via z = z{x) = {V{x) — p)/kBT. However, we 
shall usually suppress the z- and x-dependence in our notation for simplicity and just write F^{s). 

Let us see under which conditions this thermodynamic equilibrium state is stable against mechanical perturbations. 
Displacing a volume- element of fluid mechanically in the direction of increasing pressure, i.e. in the direction of /, 
its volume is compressed adiabatically, so that its density is increased, per unit displacement, by {dpo / OPq) Po, 
whereas the density in the ambient equilibrium-gas changes by Vpo in the same displacement. A restoring force per 
unit volume in the direction opposite to the displacement 



/ 



{dpQ/dPo)sVPo - Vpo 



< 



(2.6) 



must result for a n ie chan ically stable state F] U sing the relation {dpo/dPo)s = 3po/5Poi valid for the ideal quantum- 
gases, and eqs.(2^)-(2^) we may rewrite ( P^q ) as 



/ 



f2 



(2.7) 



'^A restoring force does not result if the whole fluid-layer on an equi-potential surface is displaced in the same way orthogonal 
to the equi-potential surface; because then no ambient fluid remains, which could give rise to the buoyancy force (2.6). Instead 
a new mechanical equilibrium is reached. This mechanism gives rise to the zero-frequency modes discussed later. 
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This condition is satisfied for bosons for z = {V{x) — ii)/kBT > 0, and for fermions for all positive and negative values 
of z. 

The hydrodynamic equations linearized in u, Sp, SP are then given by 



po{x)dtu = —^5P + Spf{x 
dtSP 

dtSp = -V • {poix)u) 



■ {Po{x)u) + '^Po{x)f{x)-u 



(2.8) 
(2.9) 
(2.10) 



Eliminating SP from eqs.(2.S), (2.9), then using the continuity equation ( ^.10 ) to eliminate also Sp, finally using that 
Vpo is parallel to /, one obtains immediately the closed wave-equation for the velocity field 



t o 

3 po[x) 



|v(v.ff)+v(i./)+H/(v.3) 



(2.11) 



which has been the starting point of previous works [py7 yil|J12(| . As already mentioned we do not find it the most 
convenient starting point for our present study as the boundary conditions on the velocity field at infinity, the nature 
of the function space formed by the solutions, the hermiticity or self-adjoint-ness of the wave-operator and hence 
the nature of the spectrum of eigenvalues all remain unclear, even if some particular solutions, e.g. those satisfying 
V X u = 0, ca n be constructed. To escape this impa sse w e th erefore proceed by eliminating the velocity field from 
eqs.( |2.8| - 2.10| ) by taking the time-derivative of (2.10) and ( |2.9| ) and inserting ( p.Sl ). We obtain after some calculation 



d^Sp = V^SP - fix) ■ VSp - (V • f{x))dp 



(2.12) 



3 po{x) 
5 Po{x) 
3 pq{x) 

It is convenient to define the function 



5 f^Poi^) 
3 V Pq{x) 



■VSP 



f{x)-Vdp- 



^/^4v.(^/(.) 

3 3 \Po[x) 



Sp 



and its derivative Fi-(z) with respect to its argument z in terms of which we can write 



Pq{x) UbT 



^ fPojx) 
Po{x) 



F^{{V{x)-p)/kBT) 



(2.13) 



(2.14) 



(2.15) 



= -mF' {{y{x:) - p)/kBT) 



Suppressing in the following the s ubs cript =F and also the argument {V{x) — p)/kBT of F and F' for notational 
simplicity we can rewrite eqs.( 2.12 ) - ( 2.13 ) as 



dlSp = V^SP - / • VSp - ( V 



d'iSP 



bknT 



3 TO 

5 fcjjT 

3 TO 



FV^SP - -F' 



Ff-WSp 



f) Sp 
2 

+ 3 
V- 



f-VSP 



f 



bkBT 



F 



(v./) 



Sp 



(2.16) 



(2.17) 



So far we have gained in simplicity compared to eq.( ^.ll| ) because we have now only two coupled wave-equations 
instead of three. More important, however, is the fact that it is clear physically that the density and pressure 
perturbations must go to zero in the limit of large distances from the center of a confining trap. It should be noted 
that the same cannot be said for the velocity field. Indeed, it is clear from (2^) that for |x| ^ oo where p[){x) 
the velocity field u is not necessarily bounded by the hydrodynamic equations. However, in spite of the improvement 
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of the formulation of the hnear ize d hy drodynamics we have achieved so far, the self-adjoint-ness of the wave-operator 
H defined by writing eqs.( 2.16 ), ( ^.17 ) in the form 



fSP 



H 



(2.18) 



remains to be clarified. Can a scalar product be found in which the operator H is hermitian? This is the question to 
which we turn next. 



B. Scalar product and hermiticity of the wave-operator 

In order to find a useful sc alar pro duct on the space of solutions of eqs.( 2.12 ) - ( ^.13| ) we proceed as follows. First 
we find a Lagrangian for eqs.( 2.1^ ) - ( 2.13| ), which must be a quadratic functional of Sp{x, t), SP{x, t). It can be found 
by making a gene ral ansat z and comparing the coefficient-functions of the resulting Euler-Lagrange equations with 
those in eqs.( ^.12| ) - ( 2.13 ). From the Lagrangian density we can pass to the associated 'energy density' H whose 
space-integral 



E = / (f'xn 



(2.19) 



must be conserved by time-translation invariance. We shall then define the scalar product (P1IP2) in such a way that 
the conserved 'energy' takes the form 



E = (P|HP) 

for vectors P satisfying the time-independent wave-equation HP = cj^P. 
Thus we begin with the Lagrangian 



(2.20) 



L = / (PxL 



(2.21) 



for whose density we find after some calculation 

^ = f {dtSPf + a (dtSP) (dtSp) - ^t^^idtSpf 



3 m 



Ff3 ] ( V(5P ) +E- { Sp\7SP - SPS/5p ) + ISPSp + lj{Sp) 



(2.22) 



where the coefficients a, (3, E, /, and J are defined by 

K 



P 
E 
I 
J 



(2 + 5F')i^T 



m l + F' 



-a 



keT F 
ia(5F' + 2)/ 

^ I a + - 



3 m 



|(5F' + 2)/^ 



The coefhcient K in the relation for a is arbitrary, and can be used for normalization. 
Next we pass to the associated density Ti. defined by the Legendre transformation 



(2.23) 

(2.24) 
(2.25) 
(2.26) 
(2.27) 



n = (dtdp) 



dC 



(dtSP) 



dC 

d{dt6P) 



(2.28) 
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with the conserved quantity (2.19). It is now useful to employ the vector-notation already defined in eq.(2.18) by 
defining 



Fix,t) 



6P{x, t) 
Sp{x,t) 



With the harmonic time-dependence 

P{x, t) = P(f) COs{ijjt + Lp) 

the conserved quantity E can be written as 



(2.29) 



(2.30) 



E = 0}"^ fiiv?{ut + ip) j 



+ cof?{ujt + (p) d^x 



Pfi:r,\2 , rnr 5 ksT Fa 2 

-iSP)^ + a5PSp---- — iSpr 



a 



3 m 



(2.31) 



E ■ {5pV5P - 6P^Sp^ - UPSp - ^Ji6p)^ 



In order to meet our g oal (2.2C) we should define the scalar product (P1IP2) in such a way that the coefficients of 
sin^ and cos^ in (2.31) become both equal to Ci;^(P|P). Therefore, we can conclude from the coefficient of sin^ that 
the norm becomes 



(P|P> 



d^x 



^5P^ ^ a5P5p-l^Fa{5pf 
2 6 m 



(2.32) 



Using the relations ( 2.23| ), (2.24) for a and (3 it can be checked that the norm is positive, as required, if a < 0, (which 
can always be achieved by the choice of the constant K), and /3 > 0, which requires the inequality 1 -|- F' > 0, and 



-|-^^i^a/3 — > which in turn requires the stronger inequality 



bF'{{V{x)- p)/kBT) + 2> 



(2.33) 



Using the definition ( 2.14 ) of = F^ it is easy to check that (|2.33 ) is equivalent to the stability condition ( ^.6|) of 
the thern iodyn amic equilibriimr state. The functions F(z), F'(z) are plotted in figs.|l| for bosons and |^ for fermions. 
From eq.(2.32) we deduce that the scalar product on the complexified space of solutions can be defined as 



(P1IP2 



J d^x{5P;,Spl)S 



SP2 

Sp2 



with the matrix 



S= ' 2 2 



a _5 ksT Fa 
2 3 m 2 



The coefficient of cos in eq.(2.31) can now be checked to be of the form 



(P|HP) 



(PIP) 



(2.34) 



(2.35) 



(2.36) 



for vectors |P) satisfying H|P) = llj^\P). Indeed, restricting to a real space, for simplicity, because H is real, we find 
by direct evaluation 



(P|HP)= J d^x{6P,6p)S-H(^^^^ 



= / d'^x 



- a + -- 



3 m 



F/3 j {v6Py - E ■ {SpVSP - SPVSp^ - ISP6p - ^{Sp) 



(2.37) 



We can furthermore show by direct calculation that the wave-operator is hermitian in the scalar product (2.34), if 
suitable boundary conditions are imposed. We find after straightforward, but lengthy calculation 
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(PilHPa) - (P2|HPi>= J rf^x (PiS • HP2 - P2S • HPi 



(2.38) 



i6PiWSP2 - SP2VSP1) - fiSPiSp2 - SP2Spi)\ } 



Thus, we must impose boundary conditions at infinity in such a way that the surface-integral we obtain from ( 2.3 



by the application of Gauss' law vanishes. Since the coefhcient functions a and (3 grow for \x\ 00 like exp ' — tL 



the fluctuations 6P(x) and 6p{x) for all solutions must vanish sufficiently rapidly for \x\ — > 00. 
Finally, let us transform the scalar product (2.34) to the more symmetrical form 



(P1IP2 



d^x {ulu2 + W1W2) 



by the linear transformation 



with a matrix M in lower triangular form 



fSP 



M 



M 



aO 
b c 



(2.39) 



(2.40) 



(2.41) 



diagonalizing and normalizing the kinetic term in the Lagrangian C (and hence also in the 'energy-density' Ti). The 
latter requirement yields 



m y a{5F' + 2) 
9m 1 



SfcsT ^-aF{5F' + 2) 
3rn 1 



In the new variables u, v the wave-equation now reads 



(2.42) 



u 



with the manifestly hermitian wave-operator 



H 



/ 5 kBTV-FV 
3 m 



(M 

6 



eOksT F ■ 



V5F'+2f-V I N/3m V5F'+2p 
V3 SOksT F ■ 



VV5F'+2/ 



V3 



73m V5F'+2p \ 
SOksT F ■ 



m {5F'+2)p 
bksT F 



(2.43) 



(2.44) 



C. Zero- frequency modes and isothermal modes 



For arbitrary temperature T and trap-potential V{x) eqs.(2.12), (2.13) possess a class of exact time-independent 
solutions, which depend on an arbitrary function G{V{x)), and its derivative G" = dG/dV , namely 



5p{x) = -eG'{V{x)) 

SP(x) = —G(V(x)). 
m 



(2.45) 
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e is a parameter which is sufficiently smah to make the hnearized theory consistent. The norm (2.32) of these solutions 
|Po) is 



(PolPo) 



2m 



1 + F' 
ksTF 



+ 2GG" + -kBTFC'"^ 



(2.46) 



Since all functions under the integral depend on x only via V{x) the integration J cPx... can be replaced by 
const J dVVV..., if V is a homogeneous fun ction of second order of x, e.g. a parabolic potential. The scalar product 
exists and is positive under the condition ( 2.33 ), if G{V) vanishes sufficiently rapidly, e.g. like exp(— l^/fc^T), for 
V ^ oo. Then these solutions belong to the Hilbert-space and have to be considered. Physically, they appear because 
of the coexistence of a continuum of mechanical equilibrium states and the thermodynamic equilibrium. As we shall 
see in the following section these states are not isolated from all the other states but occur, for any local wave-number, 
as the end-point of a spectral branch of states if the local wave-n umbe r is tu rned in the direction of / = — ^VV. 

There are some further exact solutions of the wave-equations (2.16, 2.17) which hold for all temper ature s in the 
fermionic -and all temperature T > Tc in the bosonic case. They are obtained by extending the ansatz (2.45) for the 
zero-frequency modes according to 



fSPix,t) 
\dp{x, t) 



;^G(F(x)) 
~G'(y{x)) 



(2.47) 



with a, /3, 7 = or 1. Inserting this ansatz in eq.( p.l6 ), and using the property that for a = (3 
zero-frequency mode, we find after a simple calculation 

(w^ - aul - I3ujI - -iujl)G'{V) = 0. 



7 = eq.(2.47) is a 



(2.48) 



Next we insert the ansatz also in eq.(2.17) and obtain by a similar calculation 

(^u^ + {au;j + Pu:l + ^u:l){^-F' + ^)^ G{V) + ^fcsT(ac.? + /3c.| + ^u:l)FG'{V) 
Eq.( ^.48" ) determines the mode- frequencies as 



= 0. 



(2.49) 



U)af3-f 



(2.50) 



while ( 2.49 ), for a,/3,j not all equal to zero, fixes the yet undetermined function G{V) in the ansatz ( 2.47 ) as 



GiV)^constFil,^ 



-). 



(2.51) 



It follows with ( 2.47 ) that SP{x,t) and 5p{x,t) for these modes are related by 

SP{x, t) 



1 rr Xl^hji] 



V{S) 



2' feflT 



Sp{x,t), 



(2.52) 



which is the relation between changes of pressure and density implied by the local thermodynamic equilibrium ( |2.2| ) 
if the temperature is kept constant. These isothermal modes were already found in |^ for the special case of isotropic 
and axially sy mmet ric parabolic trap-potentials. 

The modes ( 2.50 ) contain as special cases the three Kohn-modes ujiqo — cji,Li;oio = '.^2,'.^ooi = ^^3, corresponding 
to oscillations of the center of mass of the trapped gas. It is interesting to note that coUisionless Kohn-modes of the 
form 



(2.53) 



with the same frequencies Ui also exist. It therefore follows from the phenomenological formula (1.7) that these modes 
are not damped by the relaxation mechanisms present in the system, in agreement with the general statement made 
by the Kohn-thcorcm. 

For fermions the result ( ^.5C| ) and (2.51), (2.47) for the frequencies and mode-functions apply to all t empe ratures 
and can therefore also be extrapolated to T ^ 0. Indeed, for T ^ 0, modes with the frequencies (2.5C) where 
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found in 1^. In order to compare also the mode- functions we use the Bethe-Sommerfeld expansion to evaluate the 
Fermi-integrals asymptotically -F+(s, z) {—zy/s\ for z — > — oo and find for T — s- 



3 m 



(2.54) 



in agreement with 

Finally, another set of solutions for arbitrary temperature is found by generalizing an ansatz of Q for a trap without 
axial symmetry by putting 

u,{x,t) = AiX.e-''^* (2.55) 
with three constants Ai. Inserting directly into eqs.(2.9, 2.1(]| ) we obtain the density and pressure modes 



im f mkBT 
dp{x,t)= ~— -— 2- 



sp{x, ty 



ikfiT ( mknT 



3/2 



3/2 



2/ knT 



,2 2 



2 2 



(2.56) 
(2.57) 



A comparison of (2.56)-(2.57) with the local equilibrium relations (|l.4[)-(1.5) reveals that the temperature oscillates 



in these modes with a spatially constant amplitude proportional to Finally, using the results (2. 56), (2. 57) 

in the equation (2.S) for momentum conservation we arrive at the eigenvalue problem: 



(2.58) 



The eigenvector A and the eigenvalue ijP are clearly temperature independent and follow from the cubic secular 
equation 



20 



(2.59) 



In the special case of an axially symmetric trap the cubic equation can be reduced to a quadratic one and a result 
first obtained in |^ is recovered. 



III. SHORT WAVE-LENGTH SOLUTIONS 



The two coupled wave-equations derived in the previous section in various forms are difficult to solve for arbitrary 
temperature in a system which is made spatially inhomogeneous by an external potential V{x) 7^ 0. An exception, 
however, are waves of wave-lengths, which are short on the spatial scale on which V{x) and hence also Po(S')j Po(2^) 
vary. Such waves, in the representation with u(x, t), v{x, t), can be written as 



u{x, t) 
v(x, t) 



ba{x) 



(3.1) 



The eikonal S{x,t) defines the local wave- vector by the relation 



k{x,t) = V5(f,t) 



(3.2) 



The amplitudes ao,6o and also k vary slowly in space, on the same scale as V{x), Po{x), p(){x). The frequency uj is 



independent of x, t. Inserting the ansatz in the equation (2.43) and neglecting derivatives of oq, bo and fc, and assuming 
also that ^ |V • /|, we obtain the secular equation as the vanishing of the determinant 



5 fcflT 
3 m 



pp , ±_m_Il _ 2 



^VW^{k ■ f) 



V5F'+2 m A 



30 



-^^5F' + 2{k- f) 



n/3 V5F'+2 
30 F 



3I J 



?2 



5F'+2 m f2 _ 2 
5F koTJ ^ 



(3.3) 
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We note that terms with and k ■ f are essential to keep in this approximation together with the fc^-terms, because 
I/I grows at large distances from the trap-center, at least for parabolic traps, and provides the physically crucial 
confining mechanism. On the other hand, ao{x),bo{x), and V • / do not grow in a similar way and are therefore 
consistently neg ligible . 

From eq.(3.3) we deduce the local dispersion-law for waves of short-wavelength 27r/fc 



— uJ±{k, x) 



Fk^ 



kgT \60F 



1 



5F' + 2 



f 



bknT 



Fk^ 



knT V 60F 



1 



hF' 



n 2 



--(5F' + 2) (/xfc^ 



In the same level of approximation the pressure and density oscillations are related by 



5P{k,x) 



1 r 



if-k 



5p{k, x) 



(3.4) 



(3.5) 



as follows from eg. ( 2. 12 ). 

The dispersion law ( p.4D contains a lot of physics and will be discussed now. First we note that for f ^ there are 
two branches of the dispersion law, one of high frequency and another one of lower frequency, which are both physical. 
Thus, there are two different types of waves in these systems. Both branches correspond to frequencies ut^ > for all 
k, i.e. to stable oscillation waves. Another simple observation is that the local dispersion-relation is anisotropic and 
depends on the angle between / and k. The physical nature of the two branches is most easily seen by assuming that 



the angle between / and k is sufficiently small to permit the expansion of the square-root in (3.4) in the second term 
of its radicand. We obtain then to lowest non-vanishing order 



Fk^ 



( 1 



5F' 



keT V60F 



f 



(3.6) 



5F' 



(/X ky 



3 5 fcflT pi 2 I rn ( 1 . 5F'+2 \ f2 
3 m ^ ^ ksT KeOF 5 )J 



(3.7) 



The high-frequency branch is easily recognized in this limit as an adiabatic sound mode, in particular if the identity 
for the adiabatic sound velocity 



dPo 
dpo 



5 Pn 5 knT . 



3 Po 



-F 



(3.8) 



is used, which is valid for the ideal quantum gases, with the definition of F = Fip by eqs.(2.14), (2.15). The low- 



frequency branch reaches its lowest frequency uj- = 0, in the present approximation, for all waves traveling locally in 
the direction parallel to the force of the trap /, so that / x fc = 0. The existence of such zero-frequency modes has 
already been discussed in the preceding section. Looking at exact solutions at high-temperature in the next section 
we shall see them appear again. 

The low- frequency branch for given |fc| achieves its highest frequency if the wave propagates locally in directions 
orthogonal to /. The maximum frequency for modes orthogonal to / is then reached for short wave-lengths 



/ 1 



keT V 60F 



/ 



and given by 



5^^' + 2|/| 



(3.9) 



(3.10) 



Waves with the properties of the low-frequency branch solutions found here are typical for media which are stratified 
by an external force and are called 'internal waves'. One of their surprising and counter- intuitive properties is that in 
regions where ( ^^ ) applies the group- velocity VfeW_(fc) is orthogonal to the wave- vector k. For a text-book discussion 
of such waves see [|l^. Indeed, the dispersion-relations (3. 6), (3. 7) can e.g. be directly compared with eqs.(53),(54) 
given there. 
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IV. SOLUTION IN THE HIGH-TEMPERATURE REGION 



A. Hilbert-space of polynomial solutions 



The high-temperature regime is defined by 



deg 



(4.1) 



where Tdeg is the degeneracy temperature at which the de Broghe wave-length becomes of the order of the mean 
particle distance. In this regime we can approximate 



(4.2) 



and poix) = pa{0)e v(^)/^bT ^ Po{x) = ^^po{x). Is is useful to introduce new dimensionless variables Pi{x), pi{x) 
via the definitions 



SP{x,t) - -^po(x)Pi(f,t) 

TO 

Sp{x,t) = pa{x)pi{x,t) 
and to rewrite the coupled wave-equations for 5P, Sp as coupled equations for Pi and 

ST 



Pi-pi=Ti = 



T 



(4.3) 



(4.4) 



where 5T[x,t) is the deviation of the temperature from equilibrium. Separating the time-dependence e we arrive 



at 



Pi 



( ksT 5y72 
m 3 ^ 



/•V, -![(v-/)+/-vl-^/^\ 



V 



rn 3 ' 



,Pi - pi, 

We specialize these equations to a harmonic potential 



V-/ +/-V 



Pi 

Pi -Pi 



(4.5) 



T//' ^\ '"■ / 2 2 I 2 2 I 2 2\ 

V{x)^ — {uj^x +ujyy +uj^z) 



(4.6) 



It is then clear, via term by term inspection, that there are polynomial solutions of eq.(45), in which Pi is polynomial 
in the Cartesian components of x of total order n and Pi — Pi is a corresponding polynomial of total order n — 2. This 
is so because then each term on the right-hand side of (^|^) either decreases the total order of the polynomial on the 
left-hand side by 2 via the operation of or keeps the same order of the polynomial. It follows that we can pick 
freely as the highest total power an integer n for Pi and determine the eigenvalues by comparing the coefficients 
of all terms with this highest total power, imposing the condition for nontrivial solvability. It is clear that we ge t 
polynomials of arbitrary total order in this way. Moreover these solutions lie in the Hilbert-space because via eq. ( [1.3| ) 
5P and 5p fall off sufficiently rapidly for |a;| — > oo, if Pi and pi are polynomials in the Cartesian components of x. 

To see this explicitly we specialize the scalar product (2.34) for the present high-temperature case by replacing 
F{i/2,{V{x)-p)/kBT) = exp{-{V{x)-p)/kBT), = 1, F!^ = which gives a = -{K/2) exp{{V{x)- p)/kBT), f] = 
— (m/kBT)a, and the scalar product 



(P|P) = ?¥^(Po(0))^- 



d-^xe 



Pi Pi - Pi pi 



(4.7) 



It is also clear now, that for harmonic potentials V{x) the polynomial solutions for Pi and pi of all orders are complete 
in a space with this scalar product, a fact which is very familiar from the quantum mechanics of the harmonic oscillator. 

To be specific let us consider a polynomial of x, y, z for Pi of total order n. It has then terms of highest total order 
of the form 



Pi 



E 



lower order 
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with ni+n2+n3 = n, and ("^^) = (n + 2){n + l)/2 different coefficientsA„j„2„3. For the same mode Pi — pi must be 
a polynomial of total order n — 2 with (2) = ^ l)/2 terms x'"-^y"^z"^ of highest total order ni + n2 + n-s ~ n — 2 



with coefficients i?nin2n3- 



The solvability-condition for the linear homogeneous equations connecting all these coefficients then gives a secular 
equation for the eigenvalues of order 



n + 2\ (n 



+ 1 (4.9) 



with just as many solutions, (n + + 2)/2 of these modes can be considered as modes of sound waves, modified 
by the external potential, while the remaining n(n — l)/2 can be considered as modes of internal waves, modified by 
their coupling to sound waves. The second kind of modes therefore exists only for n>2. 

How many linearly independent zero-frequency modes appear at a given total order nl If n is odd, the polyn omial s 



for P\ and p\ have odd parity and cannot describe a zero- frequency mode, which must have even parity by eq.( 2.45 ). 
If n is even, on the other hand, there is precisely one linearly independent zero-frequency mode associated with that 
integer, which may be written in the form 



n/2 



The order of the secular equation and the number of frequencies with given value of n grows quadratically with n. 



For n we have just one coefficient Aqoo and P\ — p\~ 0, and the frequency uj^ = 0. We shall see in section IVD 
that this mode must be excluded because of particle-number conservation. Thus, there is actually no physical mode 
with n — Q. For n = 1 we have already three coefficients Aioo, ^oiOi ^001 but no B-coefficient yet. The equations for 
the three A-coefficients are decoupled and we obtain the three frequencies 

^^100 = t^a:, WoiO = t^a, ^001 = (4-11) 

corresponding to a rigid center of mass motion of the atom-cloud in the trap. These are the high-temperature limits 



of the Kohn-modes already encountered in section II C . Also the other temperature-independent modes discussed 
there appear again as follows: For n = 2 we have 6 A-coefficients and one B-coefficient and all together 7 solutions 
for Because of separate symmetry of the trapping potential under reflections of the x-, y-, and z-axis the three 
equations for Ano, Aioi, Aqu are decoupled from each other and the rest and give immediately 

'^iio = ^a;2 -hw^^wioi = _|_ w2,u;oii = ^ujI + ujI- (4.12) 
Similarly, for n — Z the coefficient Am is decoupled from the rest, and the corresponding mode frequency is 



ujm^ ^iol+^l+ujl. (4.13) 
All these decoupled modes are isothermal modes, corresponding to 5T = 0. In fact, we have already seen in subsection 



[IC that these modes exist for all temperatures (for bosons for T > Tc). The remaining four coefficients for n ~ 

2, ^200: ^0207 ^002 and i?ooo are coupled. Of the four resulting coupled modes one is a zero-frequency mode lo^^ — 0, 
for which 

1 4(0) _ 1 4(0) _ 1 4(0) 
-^200 - ^^20 - ^^02 - ^^000. (4.14) 

This is the only internal- wave mode with n — 2. The remaining three modes with n — 2 are modified sound modes, for 



whose squared freque ncies uJ^ the cubic equation. (2.59) is re-obtained as it must be, since the corresponding modes 
were found in section II C to exist for all temperatures with temperature-independent mode-frequencies (but with 
mode- functions for pressure and density which depend on temperature and the quantum-statistics). The solutions 
of the cubic secular equation are simple if one trap frequency, say ojz, is much smaller or much larger than the two 
others. Then one solution for liP is small, 

LO-d—LOz, (4.15) 

or large 



12 



respectively, while the other two are large and given by 



or small 



W4 



± 



2^2 



y ' ^ xf 



15 



96 



2,, ,2 



(4.17) 



(4.18) 



respectively. Other simple solutions are obtaine d in t he isotropic case and the axially symmetric case, cf. the following 
sections. Quartic equations for cj^ similar to eq.( 2.59| ) arc obtained from the solvability conditions for the three quadru- 
ples of amplitudes of the n = 3 modes (A102, ^120, ^300, -Bioo), (^012, ^210, ^030, -Boio), and (A021, ^201, ^oos^ooi), 
which decouple from each other because they differ in parity. The solutions of these equations and of (2.59) for 
arbitrary trap-frequencies are tedious expressions and will not be given here. They and the eigenvalues for n > 3 
can easily be determined numerically for specific ratios of the trap-frequencies when needed, where efficient use can 
be made of the already mentioned conservation of the x-, y-, and z-reflection parities. 



The solution for the isotropic case lox 



B. Isotropic case 

cjo in the high-temperature regime has already been given by 



Bruun and Clark ||7|] based on the velocity equation (2.11). It is not clear, however, how the Hilbert space is defined in 
terms of the velocity variables. Indeed, the velocity-field turns out to be polynomial like p\ and P\ and is therefore not 
a square-integrable field. It is therefore worthwhile to check this case again, using our present description. Introducing 
dimensionless frequency and space-coordinate via 



(4.19) 



and imposing a polynomial ansatz 

Pi(r) =r'Fz„,(e,(p)Q„(r2) 
Pi(f) - pi(f) = r'li„,(e, (p)r„_i(r2), 

with Qn and T„ polynomials of order n and y;,„(0, Lp) the spherical harmonics, we obtain 



2 / 

' 3 



2(; + i) d 

r dr 



2 



-2- 



21 2d 
3" ^ Vd^ 



(4.20) 

(4.21) 
(4.22) 



We should mention that the quantum- number n defined by the polynomial order of Q„, T„ differs from the quantum- 



number n defined in section IV A, which corresponds to 2n + I within the present definitions. As was mentioned 



already, the characteristic equation follows from the linear set of equations for the largest powers. Let us consider 
separately 
(i) n = : 



Eq. ([4.22D is consistent with T_ 



and gives no additional information. From eq.(4.21) we conclude 

(f]2 _ /)4") ^ 0. 



(4.23) 



(4.24) 



Thus we have the spectrum 



13 



UJ = VIluq 



(4.25) 



for the modes 

mm 

ST = 0. 

These isothermal modes were first found by Griffin, Wu and Stringari |^ 
(ii) n > 1: 



(4.26) 



From (4.21) we obtain for 



z ^ 0, n : 
i = 0: 



From (4.22) we get correspondingly 

4 



i=Q 
n-1 



i=0 



in' -I- 2n)4") + bIc\ = 
{n^-i- 2z)a|") + B^^, - ^n^Bl"^ = 



i{2i + 21 + 1)A["' + 



2 2(Z + 1) Ai 



= 



(4.27) 

(4.28) 
(4.29) 

(4.30) 
(4.31) 



Let us consider the highest power i = n first, which yields two linear homogeneous algebraic equations for B^2i 
of the form 



(4.32) 



|n [2n + 21 + 1] n^- - ^ 

-l-2n 1 



The solvability condition gives for n > 1 



n^--n^{l + 2n+-) + -l{l + l) = (4.34) 
3 5 3 



"1,2 



■'1. 2 




? + 2n + 



-Ul + l) 



(4.35) 



This result was first obtained by Bruun and Clark 0. We differ from their result only for the case n = 0, where we 
f ind t hat only one branch of the solution ( 4.35| ) exist s, nam ely the branch which yields — I, as seen from (4.23), 
(4.24). In the case n > 1,Z = one of the solutions ( 4.35| ) vanishes. It gives the zero-frequency solutions discussed 



in previous sections. The other solution at fl'^ — w^/wg = |(2n + |) behaves normally and describes a sound-mode, 
as can e.g. be seen from the factor 5/3 which is characteristic of the inverse of the adiabatic compressibility in ideal 
gases. 

Let us turn to the mode- functions in the case n > 1. Their coefficients can be obtained recursively starting with 



Bq^' , which is fixed only by normalization, then solving (4.30) 
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(4.36) 



then rising to higher values of i by first solving for 

1 sn^ - 2{l 



A 



in) 



-B 



(n) 



4 i(2i + 2/ + l) 

and then using (4.29) with (4.37) and ( f.34| ) to obtain also _B,|"'' in terms of as 



B. 



(") 



i{2i + 2l + 1) 



-B 



(n) 



Solving the recursion-relation (4.3^) we obtain 



<n) {l-n){l-n + l)---{l-n + i- 



1) I.SH 



(21±3) . + 1) . . . (2!±3 + , _ 1) ,!2^ 

which identifies the polynomial T„_i(r^) as the confluent hypergeometric function 



r„_i(r2) ^ ^ V^^ = <\Fi(1 - n, 

-(i+i) 



j=0 



2Z + 3 r2 



(4.37) 



(4.38) 



(4.39) 



(4.40) 



which is proportional to the generalized Laguerre-polynomial i„_i (^)- Substituting the result into eq.(4.37) we 
obtain also the polynomial Qn{r^) as the combination of hypergeometric functions 



g„(r2) = ^A 



B, 



(«), 



i=0 







21' 




\- 






li^i ( 


[An 


' 3" 





i^^i -n, 



-n, ■ 



2/ + 1 



2/ + 3 



21 + 1 
2n~ 

which is proportional to (Snl ^ + 2l)L^n'^^\'-^) - 2(1 + 21 + 2n)Lll^^\^). The physical modes are then 



SP{x, t) 

5T{x, t) 
T 



Po(x)r-'ri„(e,<^)Q„(r2) 



r'r,,„(e,(^)T„_i(r2) 



(4.41) 

(4.42) 
(4.43) 



It is quite remarkable that the spatial perturbation of the temperature in the two physically very different branches 
of the spectrum, the sound- modes and the internal modes, is exactly the same, because T„_i(r^) is independent of 
and therefore the same for both branches. (We note, however, that this is strictly true only in the Boltzmann 
limit). The spatial distribution of the pressure and the density, on the other hand, is very different for both kinds of 
modes as one would expect. 

For I = and arbitrary n we obtain the mode-functions of the zero-frequency modes. They form the bottom of a 
ladder of rotational modes for each value of n. The mode-functions are given by combinations of Hermite-polynomials 
of the radial variable {mojQX^ / UbTY^'^ and form a complete set in the Hilbert-space of radial functions defined by the 
scalar product (4.7). This simply means that an arbitrary radial, i.e. angle-independent, mode-function within our 
Hilbert-space is a zero-frequency mode. 



C. Anisotropic harmonic potential and conserved operator 



Let us now turn to the case of an anisotropic harmonic trapping potential. Then the generator of rotations 
IL = —il{x X V) no l onger commutes with the wave-operator H defined by the matrix-differential operator on the 
right-hand side of (4.5). On the other hand, for harmonic trapping potentials the polynomial solutions discussed in 
section p^V A continue to exist also in this case. Therefore, one must strongly suspect that a complete set of operators 
commuting with H exists also in the fully anisotropic case. For Bose-Einstcin condensed bosons and for fermions at 
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temperature T = a similar situation prevails (but not at temperatures between the high- and the low-temperature 
limits, cf. the discussion below), and two operators commuting with the corresponding wave-operator were constructed 
in our previous papers [Q. One way to find these operators is to introduce elliptic coordinates and to separate 
the wave-equation in these coordinates. The separation-constants introduced by this procedure appear naturally 
as eigenvalues of certain differential operators, any combination of which can be identified with the searched for 
commuting operators. In this procedure was fully carried through in the low-temperature limit. 

For our present high-temperature regime we shall not present the most general anisotropic case. Rather we restrict 
ourselves to axially symmetric traps with uj^ = = "^.Lj where one of the conserved operators, namely \Lz with 

i^ = -*f-i&-2/i&V (4-44) 



dy dx 

still follows from symmetry and we only need to find a second one to demonstrate integrability. 

One can find a conserving operator via the separation of the wave equation for Pi in elliptical cylindric coordinates. 
This procedure leads to an operator B', which has the important property 

(4.45) 



where the eigenvalue B is the separation constant. B' has the form 



B'-( ^ J ^\], (4.46) 



R 



5 



where (with u'^ = uj. 



R-x W- ^ (—^ + —^ + —^] (4 47) 



Returning back to our original variables SP and Sp one can prove via explicit calculation that the operator B, which 
corresponds to B in that representation, has vanishing commutator with H: 

[H,B]=0 (4.48) 



Let us determine the spectrum of B in the Hilbert-space of polynomials discussed in section IV A , whe re Pi is a 
polynomial of total order n and Pi — pi a polynomial of total order n — 2. We obtain directly from (4.45), (4.47) the 
eigenvalue 

4 

B = n- - (4.49) 
5 



i.e. B is the conserved operator which introduces the principal quantum number n of section [V A which is the total 
number of nodal surfaces of the pressure and density oscillations in any given mode. 



D. Restrictions by particle- number conservation 



In the Hilbert-space defined by the scalar product (2.32) there may still be unphysical states, which do not satisfy 
the conservation of the total number of particles 



N : 



1 



(fx {po{x) + Sp{x, tj) 



(4.50) 



which applies to a closed system. In terms of the variable pi — Sp/po eq.(4.50), with N = (fixpo{x) and 
Po{x) = po(0)e~^^^^/'^^"^, implies the condition 



d-^xe 



^V(x)/kBT 



Pi =0 



(4.51) 



In order to see which modes satisfy (4.51) we use the existence of the conserved B with ( [1.46 ), from which, by 
elimination of Pi , the eigenvalue equation for pi 
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R PI 



can be derived very easily. We can use this to insert R in front of pi in eq. (|4.51| ) to obtain 



-V{x)/kBT 



Rpi 



B 



-V(x)/kBT 



Pi 



(4.52) 



(4.53) 



Now we use the exphcit form of R and the quadratic form of the potential on the left hand side of ( 4.53| ) and apply 
partial integration to let the derivatives act on the exponential factor. Boundary terms are not incurred by this 
operati on be cause pi is polynomial and the exponential factor vanishes rapidly at infinity. We find that the left hand 
side of ( 4.53 ) vanishes, and therefore 



4 

^+5 



(4.54) 



Hence, all modes for which the eigenvalue B satisfies 



(4.55) 



satisf y the restriction ( 4.51 ) imposed by particle- number conservation. Looking at the spectrum ( 4.49 ) of B we see 
that (4.55) is satisfied for all modes wit h the exception of the mode with n = 0, which is a special zero-frequency 
mode. For this mode pi{x) = const and (4.51) is obviously violated. This mode must therefore be excluded from the 
physical spectrum. All the other modes, and in particular all the other zero-frequency modes, satisfy particle-number 
conservation. 



E. Solutions for axially symmetric traps 

We introduce scaled cylinder coordinates r,z,ip via 



knT 



Xi 



Muji 



■^r cos if, 



X2 = 



MloI 



r smip, 



X3 



(4.56) 



and the ansatz 



Pi~pi= rl'"le*™'^z"T(p2,z2) 



(4.57) 



where a = 0, 1 determines the parity under inversion of the z-axis and m = — cxd, • • • , 0, 1, 2, • • • c» is the quantum 
number of angular momentum around the z-axis. The equations to be solved are then, with f2^ = w^/w^, A = w^/w^ 



n^T = — 



2 

+ 3 



92 92 2|to| + 1 d 2a d 
dp^ dz^ p dp z dz 

d d 

\m\ + 2 + X{a + 1) + p^ + \z — 
op oz 



Q 

T 



(4.58) 



2 



T 



11 ^ d B 

TO + Aa -I- p- — h \z — 
op Oz 



Q. 



For Q and T we make a polynomial ansatz in p^, of order n and n — 1, respectively. The eigenvalue of B' for these 
solutions isi3 = 2n-|-|TO|-|-Q;— |. The terms of Q and T of highest order read 



i=0 
ri-l 



2;^2(«-i-/) 



lower order 



lower order 



(4.59) 



i=0 
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Inserting the ansatz and comparing coefficients of the highest order terms we obtain the two sets of equations 



n^Ao= [\m\ + A(2?i + a)] Aq - X^Bq 

n^Ai = -Bo + [\m\ + 2 + X{2n - 2 + a)] - X^Bi 



(4.60) 



n^Ai=-Bi^ 



21 + A(2n -21 + a)] Ai - X^Bi 



VL^An= -Bn-i + [\m\ + 2n + Xa] A^ 



and 



{n-l){2n-2l - 1 + 2a) Ai 

[\m\ + 21 + 2 + X{2n - 2^ + a - 1)] Bi 



(4.61) 



--(/ + l)(|m|+; + l)^,+i 



Let us examine simple special cases. 

For n = we have Q = Aq,T = {) and obtain 



\m\ + Xa (|m| = 0, 1, • • • ;a = 0, 1) 



(4.62) 



This result is a special case of the result ( 4.11 ) (for m — Q^a — \ and m — ±1, a ~ 0), and ( 4.12| ) (for m = ±1, a — 1 
and one mode with \m\ = 2,a = 0), and ( 4.13| ) (for one mode with \m\ — 2,a — 1), but the case \m\ > 3 has no 
simple counterpart in the fully anisotropic case. 



For rt = 1 wc get already the three coupled equations 



X^Bo 



- ^{\m\ + 2 + X + aX) 



Bo 
Bo 



[1^2 - \m\ - {2 + a)X]Ao = 
{n'^ - \m\ -2- Xa)Ai = 

^{l + 2a)Ao + ^i\m\ + l)Ai 



Specializing further to m = and a = we find the exact solutions 



with 



Ao = XAi = -Bo 



and 



(4.63) 



(4.64) 



4A + 5 ± v/l6(A- 1)2 +9 



These mode- frequencies were first obtained in 



(4.65) 



The zero- frequency mode is the a xially symmetric counterpart of ( 4.14) while (4.65) corresponds to two of the 
three mode- frequencies which solve ( 2.5£ ). The third solution of ( 2.59 ) in the axially symmetric case is uj'^ — 2(jj^ 
and is in the completely anisotropic case the counterpart of the second of the two modes described by (4.62) with 
\m\ = 2,a = 0. 

In order to account for all modes in our comparison between the axially symmetric and the f ully a nisotropic 
case, we may finally note that the second of the two modes with |m| = 2, a = 1 described by (4.62) has as a 
completely anisotropic counterpart a particular solution of the quartic secular-equation for uj'^ for the amplitudes 
(Ao2i, ^201, ^003, -Booi) appearing in the scheme of section IV A for n — 3. 
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V. NUMERICAL DETERMINATION OF THE TEMPERATURE-DEPENDENT MODE-SPECTRUM 



At intermediate temperatures the coupled wave equations ( 2.16| ), (2.17) are generally not separable and the spectrum 
can only be found numerically. Best suited for numerical work are the wave equations in the manifestly hermitian 
form ( 2.44| ) with the scalar product in the simple form (2.39). 

The numerical analysis will be performed for the axially symmetric case, choosing for the anisotropy parameter 
A = Lo\luj\ = 8, partially for historical reasons as this was the geometry of the first TOP trap at JILA ||2^. The 
number of particles is chosen as = 10^. First the chemical potential is determin ed f or the given particle-number N 
as a function of temperature. This is done in the standard way, by integrating eq.( |l.4[ ) to obtain A^(/i, T) and solving 
for /i(Af, T). The results for our chosen set of parameters are displayed in figj^ both for bosons and for fermions, 
in the domain where /i < 0, to which we s hall restrict our attention, in the following. We may remark here that it 
follows from the form of the potential and ( |l.4| ) that /i is a scaling function 5 of A^, tD = [iOx^y^z and T of the 
form ^/^cS = S{kBT/nu}N^I^). 

In order to determine the spectrum, the wave-operator (2.44) is represented in the basis of the harmonic oscillator 
eigenfunctions with widths ^JkgT /mLol and ^k^T jmuj^ in axial and radial direction, respectively. The basis is cut- 
off at a finite size of order 100 both for u and v and the resulting finite-dimensional Hermitian matrix is diagonalized. 
The size of the finite basis is varied in order to control that the eigenvalues obtained are converged numerically. 
The truncation of the basis to a finite size introduces some spurious eigenmodes and eigenfrequencies, which can be 
distinguished however, and subsequently eliminated, by the fact that they don't converge but disappear and reappear 
somewhere else as the size of the basis is varied. 

Some of the results are displayed in figs.^^ which we now discuss. Fig.^ gives an overview, in the domain 94 < 
kBT/hu} < 195, of the spectrum of eigenvalues w^, for a gas of bosons (but very similar results not shown here are 
also found for fermions). For clarity only modes with azimuthal quantum number m = and even parity are shown. 
The basis used consisted of oscillator eigenfunctions of order 2n^ in z-direction and order 2np in radial direction, with 
integer rip + Uz < 10. All frequencies obtained by the diagonalization of the matrix of H in this basis (except for the 
spurious ones introduced by the truncation of the basis) are shown in the figure. As can be seen from fig.^ the chosen 
temperature domain extends from the BEC-temperature to the high-temperature region. In fact, we have checked 
that our numerical code applied to the case of a n iso tropic trap gives a similar result, which for kBT/hw — 195 
coincides with the analytically known eigenvalues ( 4.35| ) in the high-temperature limit with high precision. This fact 
offers us the possibility to assign the quantum numbers of the high-temperature domain to the whole corresponding 
temperature-dependent branch of frequencies. 

An obvious feature of fig.^ is the band-like structure of the spectrum, which is caused by the anisotropy of the trap: 
Because of the higher stiffness of the trap in axial direction for the assumed value of A = 8, nodes of the mode- function 
in axial direction are more costly in energy Jiu>, for sound-like modes, than nodes in radial direction. Therefore we 
can assign to the 'bands' the quantum number riz of nodes of sound-like modes in axial direction with riz = 10 in 
the highest 'band' (which can consist of a single mode only in the subspace we consider in fig.^) and ti^ = in the 
lowest. It should be noted that the two lowest bands are not split and form a single broad band. Within a given 
band sound-like modes differ only by the radial quantum number Up counting the number of nodes in radial direction. 
By our choice of a finite-dimensional basis we are restricted to modes with quantum numbers rip -\- Uz < 10. The 
presence of internal waves complicates the assignment of quantum numbers, because each pair of quantum numbers 
Uz , Up appears twice, once at higher frequency for a sound-like mode and once at lower frequency for an internal wave 
mode. Thus internal modes with high values of Uz give also frequencies in the low lying bands of fig.^. The number 
of eigenfrequencies in the 'bands' depends, of course, on the size of the basis. 

By these considerations we arrive at the following assignment of quantum numbers to the frequencies shown in figJJ. 
The mode with the largest frequency forms the 'band' on the top and must be a sound-mode with Uz = 10, rip ~ 0. 
The next lower band must contain two sound-like modes with Uz = 9,np = 1 and Uz = 9,np = 0. We see from the 
figure that this exhausts already the number of modes (=2) in the second 'band' from the top, which can therefore not 
contain an internal wave mode. Similarly, in the 3''*, 4*'', 5*'*, etc. 'band' from the top there are 3, 4, 5, etc. sound-like 
modes with Hz = 8,7,6, etc. and 3,4,5, etc. different values of Hp, respectively. In the figure we can follow this 
counting of different frequencies down to the third 'band' from the bottom, thereby accounting for all eigenvalues in 
these 'bands'. It follows that none of these 'bands' contains any internal wave modes. 

The lowest two bands, on the other hand, in particular the lowest one, contain many more different frequencies than 
the 10 and 11 different sound-frequencies with = l,nz = 0, respectively. These must therefore be considered as the 
internal waves. Among the internal waves are also the zero-frequency modes. For each value of n = Hz + np > 1 in our 
subspace there is precisely 1 zero-frequency mode, so the eigenvalue = is 10-fold degenerate in our subspace. By 
their special nature the zero- frequency modes have quantum- numbers Uz = Up. Internal modes of high n = 2{nz+np)^ 
for which rip differs only slightly from Uz have frequencies close to zero 
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The discussion of the spectrum we have given depends on our arbitrary restriction of the number of nodal surfaces 
to a finite and not very large number. If modes with an arbitrary number of nodal surfaces are permitted, then the 
mode-spectrum becomes dense, due to the existence of small internal wave frequencies from mode s with arbitrarily 
high quantum- numbers. This can already be seen from the analytically determined spectrum ( 4.35 ) for the isotropic 
case. 

Within the subspace of Hilbert space we consider here there are many internal mode frequencies also below the 
geometric mean trap frequency, as can be seen with more clarity in fig.^. Only a small number of sound-modes can 
occur in this regime. The mode frequencies displayed in figs.^, ^ have a surprisingly weak temperature-dependence 
throughout the range considered. For large T this can also be seen from the analytical results. For the sound- 
modes the velocity of sound increases with temperature roughly ~ s/ksT ^ but the wavelength of any given mode also 
increases with temperature ~ \fkBT due to the expansion of the size of the thermal cloud, so that both temperature- 
dependences effectively cancel. However, for internal modes the compensation between the speed of sound and the 
wave-length does not work in the same way, as can e.g. be seen from eq.( p.7] ). In fig.^ we consider a magnification 
of a part of the frequency spectrum (in this instance for the isotropic case)where the increase of some frequencies 
with temperature can be seen, which freely cross other levels (which must have therefore different quantum numbers 
n, Z), which are nearly temperature-independent and belong to sound- like modes. However, in the axially symmetric 
case avoided crossings can also be seen, as shown in fig.[7| for two internal wave modes at s mall fr equency below 
the geometric mean trap-frequency. This indicates that the conserved operator B of section [VC ceases to exist 



at intermediate temperatures. It would be futile, therefore, to look for analytical solutions of the spectrum in the 
intermediate temperature range, as the system appears to be non-integrable. 

Generally, the differences between the results for fermions and bosons in the region above the degeneracy tem- 
perature are qualitatively not very big. One difference due to quantum statistics can be seen in fig.^, where the 
temperature dependence of internal wave modes and sound modes is displayed for a Fermi gas. The frequencies of the 
internal wave modes curve downwards, those of the sound modes curve upwards. In the Bose gas case the opposite 
tendency is found as shown in figures ^ and ^. 



VI. CONCLUSIONS 



In the present paper we have given a systematic analysis of the hydrodynamic modes of quantum gases in a 
harmonic trap with general anisotropy in the collision-dominated non-dissipative limit. Provided the hydrodynamic 
limit is applicable the analysis applies for Fermi-gases at all temperatures and for Bose-gases at temperatures above 
Bose-Einstein condensation. Our results extend previous works by allowing for traps with arbitrary anisotropy, and 
by treating bosons and fermions side by side on an equal footing. In addition to analytical solutions in certain special 
cases and limits we also present numerical solutions in the whole temperature-domain covered by the theory. Our 
analysis is based on the reduction of the five conservation-laws for the densities of mass, momentum and energy to 
two coupled wave-equations for the mass-density and the pressure. We have constructed a scalar product with a 
positive L2-iiorm on the space of solutions, in which the two-component wave-operator is hermitian and, because of 
the stability of the hydrodynamic modes which we demonstrate, it is non-negative. However, a class of solutions with 
vanishing frequencies was found which is a consequence of the existence of mechanical equilibrium states in addition 
to the unique thermodynamic equilibrium which maximizes the entropy. A further class of exact solutions consisting 
of isothermal modes was identified among which are the center of mass modes required by the Kohn-theorem. These 
results generalize earlier results by Griffin et al. |5| and Bruun and Clark |7| to traps without axial symmetry. 

We studied the two coupled wave equations for pressure and density in the short- wavelength limit. Two different 
branches of solutions could be identified in this way, the high-frequency branch being associated with pressure-driven 
sound waves, the low-frequency branch with potential-driven internal waves. The explicit dispersion relation of the 
lower branch found in the short wavelength limit and its characteristic properties like the existence of a maximal 
frequency, the anisotropy of the dispersion-relation and the orthogonality of the local group velocity and the local 
wave-vector makes the identification with internal waves manifest and unambiguous. We also examined the high- 
temperature limit of the two coupled wave equations and demonstrated the existence of polynomial solutions by 
exhibiting a conserved operator whose eigenvalues fix the respective polynomial order, or, equivalently, the number 
of nodal surfaces of the solutions. We constructed the solutions for density and pressure also explicitly. 

The wave equations were finally also solved numerically. Surprisingly, the mode-spectrum found turned out to be 
quasi-continuous, which can be understood by the overlap of the spectrum of low-frequency internal waves of short 
wavelengths with the spectrum of sound waves with large wavelengths. Within a finite-dimensional subspace of the 
Hilbert space, defined by restricting from above the number of nodal surfaces, we found a band-like structure of 
the eigenfrequencies in a strongly anisotropic trap, where the bands are labeled by the number of nodal surfaces 
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orthogonal to the more strongly confined direction while the modes within a given band differ by the number of nodes 
in the other directions. The number of internal wave modes in such a finite-dimensional subspace is also limited, and 
these modes are then found primarily in the lowest lying band. We note that the restriction to a finite subspace of 
Hilbert space is also physically motivated, since typical excitation-mechanisms like the modulation of the trapping 
potential, will also excite only modes in a certain subspace with appreciable amplitude. 

The analysis we have presented is subject to some obvious limitations which we now discuss briefly. It is clear 
that only systems in the collision-dominated hydrodynamic limit have been considered here. For bosons a necessary 
requirement is therefore a large scattering length and a sufficiently large number-density to ensure a large cross-section 
for elastic collisions. For fermions this requires, besides a large number-density, the simultaneous trapping of several 
hyperfine states, in order to allow for the interaction of the different fermionic species by elastic collisions, which would 
be forbidden for a single species by the Pauli-principle. We have neglected throughout the spin- wave excitations, which 
can also occur in the latter systems, which is permitted because they decouple from the density-waves by symmetry 
as long as the external potential is the same for all components. The collision-rates in degenerate Fermi-gases are 
suppressed by a Fermi-blocking factor compared to the classical collision rates ||2l[] and scale proportional to (T/Tp)^. 
Therefore, at least in the low-temperature domain, it is necessary to use atomic species with particularly large positive 
or negative s-wave scattering-lengths. Another limitation of our analysis is the neglection of mean-field effects of the 
interaction in comparison to the pressure term. This seems to be a rather good approximation for the experimentally 
realized trapped quantum-gases, which behave like ideal quantum-gases to a good approximation. The most severe 
limitation of our calculations is certainly the neglect of dissipation. The reason for this restriction (which is discussed 
further in |^ for bosons, and in for fermions) lies not so much in the negligibility of dissipative effects for the 
physically excited modes, but in the particular purpose we set out to achieve in this paper, namely to give an account 
of the mode-spectrum in the whole temperature-domain covered by the theory. This goal cannot be achieved so far 
with the inclusion of damping effects, but remains an interesting aim for future work. It seems clear that with the 
inclusion of damping the zero-frequency modes we have found will turn into purely over-damped modes. However, as 
was exemplified for the Kohn- modes, the phenomenological theory behind eq.(^]^) permits to obtain also a result for 
the damping of some modes (with the result that it is vanishing for the Kohn-modes) whose hydrodynamic frequencies 
are only known in the absence of dissipation JlT| , p^ . The requirement is that the mode- frequencies without damping 
are also known in the collisionless limit. Then an estimate of the collision-time can be used in eq.(l/7) to obtain an 
interpolation between the collision-dominated and the collisionless regime including the damping due to the finite value 
of the collision-tim e [ll|Jl^ . In the past this estimate proved to be quite useful in the comparison of the experimental 
results for bosons ]22|], and it may be hoped that further results along such lines not only for bosons but also for the 
rapidly developing experiments on trapped Fermi-gases may be obtained in the near future. 
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FIG. 1. The bosonic function F-{z) = _F_(|, 2;)/F_(|, z) (full line) and its derivative F'_{z) (broken line) as a function of z. 
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FIG. 2. The fermionic function F+{z) = F+(|, a) (full line) and its derivative F!^_{z) (broken line) as a function of 
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FIG. 3. The dimensionlcss chemical potential /x/ftw in the range /u. < as a function of the dimonsionless temperature 
t = kBTjTvuj with oj = (uj^ui^y^^, for bosons (full curve) and for fermions (dashed curve); number of atoms N = 10^. 
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FIG. 4. Dimensionless squared hydrodynamic mode-frequencies (uj/lu)^ of a bosonic gas above the BEC-transition as function 
of the dimensionless temperature t, for the m=0 modes of even parity with up to 10 nodal surfaces; anisotropy parameter 
A — {u)z/lj±)^ = 8; number of atoms A*' — lO''; total size of basis 132. 
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FIG. 5. Squared internal mode frequencies of fig.^ below the geometric mean trap frequency as function of temperature, for 
the same parameters as in fig.kl 
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FIG. 6. Squared mode-frequencies increasing with temperature and level crossings for a bose-gas with the same parameters 
as in figM for an isotropic trap. 
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FIG. 7. Squared mode-frequencies of internal waves decreasing with temperature, and an avoided level crossing for a bose-gas 
in an axially symmetric trap with the same parameters as in fig.^ 
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FIG. 8. Squared mode-frequencies for a fermionic gas above the Fermi- Temperature as a function oit = fcsT/to for sound 
modes (upper part) and for internal wave modes (lower part). 
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